Dynamo Onset as a First-Order Transition: Lessons from a Shell Model for 

Magnetohydrodynamics 
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We carry out systematic and high-resolution studies of dynamo action in a shell model for magne- 
tohydrodynamic (MHD) turbulence over wide ranges of the magnetic Prandtl number PrM and the 
magnetic Reynolds number ReM- Our study suggests that it is natural to think of dynamo onset 
as a nonequilibrium, first-order phase transition between two different turbulent, but statistically 
steady, states. The ratio of the magnetic and kinetic energies is a convenient order parameter for 
this transition. By using this order parameter, we obtain the stability diagram (or nonequilibrium 
phase diagram) for dynamo formation in our MHD shell model in the (Pr M ,Reyi) plane. The 
dynamo boundary, which separates dynamo and no-dynamo regions, appears to have a fractal char- 
acter. We obtain hysteretic behavior of the order parameter across this boundary and suggestions 
of nucleation-type phenomena. 

PACS numbers: 47.27.Gs,47.65.+a,05.45.-a 



I. INTRODUCTION 



The elucidation of dynamo action is a problem of cen- 
tral importance in nonlinear dynamics because it has im- 
plications for a variety of physical systems. Dynamo in- 
stabilities, which amplify weak magnetic fields in a tur- 
bulent conducting fluid, are believed to be the principal 
mechanism for the generation of magnetic fields in celes- 
tial bodies and in the interstellar medium u L p> [ j, H , 
H, 0], and in liquid- metal systems [1, [t| [Tol fill Tl37[l3| 
studied in laboratories. In these situations the kinematic 
viscosity v and the magnetic diffusivity r\ can differ by 
several orders of magnitude, so the magnetic Prandtl 
number Ptyi = vjr\ can either be very small or very 
large; e.g., PrM — 10 -2 at the base of the Sun's convec- 
tion zone, PrM — 10~ 5 in the liquid-sodium system, and 
PrM — 10 14 in the interstellar medium. This Prandtl 
number is related to the Reynolds number Re — UL/u 
and the magnetic Reynolds number PeM — UL/rj that 
characterize the conducting fluid; here L and U are 
typical length and velocity scales in the flow; clearly 
PrM = ReM /Re- 

Two dissipative scales play an important role here; 
they are the Kolmogorov scale Id [~ v 3 ^ 4 at the level 
of Kolmogorov 1941 (K41) phenomenology Oil and the 
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in K41]. For large 



Prandtl numbers, i.e., PrM 3> 1, l 1 ^ *C id so the mag- 
netic field grows predominantly in the dissipation range 
of the fluid till it is strong enough to affect the dynamics 
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of the fluid through the Lorentz force. This behavior is a 
characteristic of a small-scale turbulent dynamo, in which 
dynamo action is driven by a smooth, dissipative-scale 
velocity field. In the initial stage of growth, called the 
kinematic stage of the dynamo, the magnetic field is not 
large enough to act back on the velocity field. Dynamo 
action can be obtained for values of PeM that are large 
enough to overcome Joule dissipation; and the dynamo- 
threshold value PeMb decreases as PrM increases [nl O] . 
Pr M ^ 10 -5 in liquid-metal flows 0, El G3] so they lie 
in the small-Prandtl- number region, PrM <C 1, for which 
the growth of the magnetic energy occurs initially in the 
inertial scales of fluid turbulence, because id <C i™ ; here 
the velocity field is not smooth and the local strain rate 
is not uniform: At the K41 level the turnover velocity of 
an eddy of size I is v(£) ~ i 1 ^ ', so the rate of shearing 

(v{t)/i) ~r 2 / 3 . 

Direct numerical simulations (DNS) arc playing an in- 
creasingly important role in developing an understanding 
of such dynamo action. Most DNS studies of MHD tur- 
bulence [H, O, EH H3| have been restricted, because of 
computational constraints, either to low resolutions or to 
the case PrM = 1; small-scale dynamos, with PrM 3> 1 
have also been studied via DNS [7fl. However, given the 
large range spanned by PrM in the physical settings men- 
tioned above, some recent DNS studies of the MHD equa- 
tions have started to explore the PrM dependence of dy- 
namo action: the range of PrM covered by such pure DNS 
studies [H, HH] is quite modest (1CT 2 < Pr M < 10). 
To explore the dynamo boundary in the (Pr M : , Reu) 
plane over a large range of Pvm, one recent study [161 ] 
has used a combination of numerical methods, some of 
which require small-scale models like large-eddy simula- 
tions (LES) or lagrangian-averaged MHD (LAMHD), and 
others, like a pseudospectral DNS, in which the only ap- 
proximations are the finite number of collocation points 
and the finite step used in time marching; yet another 
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DNS study [13] has introduced hyperviscosity of order 
8 to study the low PrM regime; by using this combina- 
tion of methods these studies has been able to cover the 
range 1CP 2 < Pry[ < 10 3 and to obtain the boundary 
between dynamo and no-dynamo regions but with fairly 
large error bars. 

We have carried out extensive, high-resolution, numer- 
ical studies that have been designed to explore in de- 
tail the boundary between the dynamo and no-dynamo 
regimes in the {Pr^ , ReyA plane in a shell model for 
three-dimensional MHD [23, EE SI- This shell 
model allows us to explore a much larger range of PrM 
than is possible if we use the MHD equations. Al- 
though our study uses a simple shell model, it has the 
virtue that it can explore the boundary between dynamo 
and no-dynamo regions in great detail without resort- 
ing to the modelling of small spatial scales. Shell-model 
studies of dynamo action have also been attempted in 
Refs. [13, E3; EE E3 but these have concentrated on as- 
pects of the dynamo problem that are different from those 
we consider here. 

Our study suggests that it is natural to think of the 
boundary between dynamo and no-dynamo regimes in 
the (-FV^ 1 , Rbm) plane as a first-order phase boundary 
that is the locus of first-order, nonequilibrium phase tran- 
sitions from one nonequilibrium statistical steady state 
(NESS) to another. The first NESS is a turbulent, but 
statistically steady, conducting fluid in which the mag- 
netic energy is negligibly small compared to the kinetic 
energy; the second NESS is also a statistically steady 
turbulent state but one in which the magnetic energy is 
comparable to the kinetic energy. Indeed, the ratio of 
the magnetic and fluid energies Eb/E u turns out to be a 
convenient order parameter for this nonequilibrium phase 
transition since it vanishes in the no-dynamo phase and 
assumes a finite, nonzero value in the dynamo state. The 
other, intriguing result of our study is that the bound- 
ary between these phases is very intricate and might well 
have a fractal character; this provides an appealing ex- 
planation for the large error bars in earlier attempts to 
determine this boundary [IE E3 ■ The analogy with first- 
order transitions that we have outlined above is not su- 
perficial. As in any first-order transition we find that 
our order parameter shows hysteretic behavior [3fj| as we 
scan through the dynamo boundary by changing the forc- 
ing term at a nonzero rate. We also find some evidence 
of nucleation-type phenomena: the closer we are to the 
dynamo boundary, the longer it takes for a significant 
magnetic field to nucleate and thus lead to dynamo ac- 
tion. We compare our results with earlier studies such 
as Ref. [3l|, which have suggested that dynamo action 
occurs because of a subcritical bifurcation. 

The remaining part of this paper is organised as 
follows: In Sec. [IT] we describe the shell model for 
MHD [IE EE Ell and the numerical method we employ. 
Sec. Mil is devoted to our results and Sec. IIVI contains a 
concluding discussion. 



II. MODELS AND NUMERICAL METHODS 

To study dynamo action it is natural to use the equa- 
tions of magnetohydrodynamics(MHD). In three dimen- 
sions the MHD equations are 

q — » -i 

^ + (tf.V)£ = vV 2 u- Vp+— (6- V)b + f, (1) 
db - 

— = V x (u x b) + ?/V 2 6, (2) 

where v and rj are the kinematic viscosity and the mag- 
netic diffusivity, respectively, the effective pressure p = 
p + (6 2 /87r), and p is the pressure. For low-Mach-numbcr 
flows, to which we restrict ourselves, we use the incom- 
pressibility condition V • u(x, t) — 0; and V • b(x, t) = 0. 

As we have mentioned above, a DNS of the MHD equa- 
tions poses a significant computational challenge, even 
on the most powerful computers available today, if we 
want to cover a large part of the (Pr^ 1 , Rcm) plane and 
to locate the dynamo boundary accurately. Therefore 
one study [l(| has used a combination of LES, LAMHD, 
and DNS to obtain this boundary. We employ a com- 
plementary strategy: we use a simple shell model for 
MHD [IE El, El that allows us to carry out very ex- 
tensive numerical simulations to probe the nature of the 
dynamo boundary without using LES or LAMHD. 

Shell models comprise a set of ordinary differential 
equations with nonlinear coupling terms that mimic the 
advection terms in and respect the shell-model analogs of 
the conservation laws of the parent hydrodynamic equa- 
tions in the inviscid, unforced limit [32J, [33]. For the 
case of MHD each shell n is characterized by a complex 
velocity u n and a complex magnetic field b n in a loga- 
rithmically discretized Fourier space with wave vectors 
k n ; furthermore, there is a direct coupling only between 
velocities and magnetic fields in nearest and next-nearest 
neighbor shells. The MHD shell model equations [IE El 
are 

= ~ vk 2 n u n + i[A n (u n+1 u n+2 - b n+ ib n+2 ) 
+ B n {u 

+ C„(u n _ 2 u n _ 1 -6 n _ 2 & n _ 1 )]* + /£, (3) 
= - rjk^bn + i[D n (u n +\b n+ 2 - b n+ iu n+2 ) 

&n+l &n — l u n+l) 

2«n-l)]* + fn, (4) 

where * denotes complex conjugation, 1 < n < N, 
with N the total number of shells, the wave numbers 
k n = fco2™, with fco = 2~ 4 , and f% and the forc- 
ing terms in the equations for u n and b n , respectively. 
In our studies of dynamo action, we set /„ = 0. The 
parameters A n , B ni . . . , F n , are obtained by demand- 
ing that these equations conserve all the shell-model 
analogs of the invariants of 3DMHD, in the inviscid, 
unforced case, and reduce to the well-known Gledzer- 
Ohkitani-Yamada (GOY) shell model [H [H[ for fluid 
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turbulence if b n = 0, V n. In particular, to ensure the 
conservation of shell-model analogs of the total energy 
E T = E u + E b = (l/2)£„(|u„| 2 + \b n \ 2 ), cross helic- 
ity He = (l/2)X)n( u « & n + u V ) n), and magnetic helicity 
Hm = l) n \bn\ 2 /k n , in the unforced and inviscid 

case, and to obtain the GOY-model limit for the fluid, 
we choose 



the K41 form ~ k n if we ignore intermittency correc- 
tions. We now introduce a small seed magnetic which 
is such that E b ~ 10~ 28 and then follow the temporal 
evolution of u n and b n that is given by Eqs.([3]) and (01. 

III. RESULTS 



D n = fc„/6; E n = k n ^ 1 /3; F n = -2fc n _ 2 /3. (5) 

The only adjustable parameters are the forcing terms and 
v and 7/. The ratio v/rj yields the magnetic Prandtl num- 
ber Pru- Grashof numbers yield nondimensionalizcd 
forces [341 but, for easy comparison with earlier stud- 
ies @, h m tin m m, we use tne nuid and ma s- 

netic integral-scale fluid and magnetic Reynolds num- 
bers whose shell-model analogs are, respectively, Re — 

Ur mJl/v, Where t\ = J2( U l/ k l)/ XX U n/ fc n) and U rms = 

VJ2( u n/ k n)/#o, to = 2ir/ki and Re M = Pr M Re. 
We use the following boundary conditions: 

An— i = An = B\ = Bn = C\ = C2 = ; 
Djv-i = D N = E x = E N = Fx = F 2 = 0. (6) 

We set N = 30 and use a fifth-order, Adams-Bashforth 
scheme for solving the shell-model equations, i.e., for an 
equation of the type 



dq 



-aq+f(t), 



(7) 



we use 



q(t + 8t) = e- 2aSt q(t - St) + 



1 



-2a6t 



24a 



x [55/(t) - 59/(t - 6t) 

+ 37f(t-25t)-9f(t-35t)], 



(8) 



where St is the time step. We have found that this numer- 
ical scheme works well for the integration of Eqs.f3|) and 
so long as N < 35 and Re < 10 9 . In all our calcula- 
tions we use St = 10~ 4 . Characteristic time scales include 
the time scale for diffusion r n = £q/j] and the large-eddy- 
turnover time tl — ^o/ u rms, where £q = 2n/ki is the 
box-size length scale and it rms is the root-mean-square 
velocity. 

The initial conditions we use are as follows: We first 
obtain a statistically steady state for the GOY-shell- 
model equations, which are obtained from Eq.Q by 
setting all b n = 0; the forcing terms are chosen to be 
fn = /o(l + i)5 n ,i, with f = 5.0 x 10~ 3 in all our runs, 
except ones in which we study hysteretic behavior, and 
/k = 0. We choose the GOY-model shell velocities at 

— 1/3 

time t = to be u n = k n exp(iip n ), with ip n a random 
phase distributed uniformly on the interval [0,27r). To 
make sure we have a statistically steady state we evolve 
the shell velocities u n till t = 5 x 10 s . This yields the 



UJ 
LU 



2 



x 10" 




0.5 1 1.5 2 2.5 



N if x1cr 



4 
t/x 



x 10" 



shell- model energy spectrum E u (k n 



- /k n that has 



FIG. 1: (Color online) Representative plots of the dynamo or- 
der parameter Ebj ' E u versus time t/T v , with r v the magnetic- 
diffusion time, in the dynamo region (blue, dashed curve), 
near the dynamo boundary (green, full line), and in the no- 
dynamo regime (red, full line in the inset). 

Given the numerical scheme that we have described in 
the previous Section, we obtain the time series for u n 
and b n from the MHD-shell- model equations. An analy- 
sis of these time series shows two types of nonequilibrium 
statistical steady states (NESS). We refer to the first as 
the no-dynamo state and to the second as the dynamo 
state. These states have been found in several earlier 
studies such as Refs. 0, El, El 0J, M, EH, H3- Our 
main goal is to explore in detail the phase boundary be- 
tween these two states. This can be done most easily by 
the introduction of a dynamo order parameter a natural 
candidate for which is the ratio Ef,/E u , where the fluid 
and magnetic energies are, respectively, E u = \ ^2, n \u n \ 2 
and Et = hY^n IM 2 - Representative plots of this or- 
der parameter are given as functions of time t in Fig. [T] 
The blue, dashed curve shows the evolution of Eb/E u 
in the dynamo regime; note that here the dynamo or- 
der parameter rises rapidly, fluctuates significantly for 
t/ T v ^ 2 x 10 -4 , and finally reaches a statistical steady 
state with equipartition, i.e., Eb/E u ~ 1. The red, full 
curve in the inset of Fig. Q] shows how Ej, /E u vanishes 
rapidly in the no-dynamo state. The behavior of the 
dynamo order parameter is more complicated than these 
two simple possibilities in the vicinity of the phase bound- 
ary between dynamo and no-dynamo states as shown by 
the green, full curve in Fig. [TJ E\,jE u rises much more 
slowly from zero than in the dynamo regime and then 
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FIG. 2: (Color online) Semi-log (base 10) plots of the kinetic (red full curve) and magnetic (blue dashed curve) energies 
versus time for (a) Pr M = 10 2 , (b) Pr M = 1, (c) Pr M = 10"\ (d) Pr M = 10~ 2 , (e) Pr M = 10~ 3 , (f) Pr M = 1CT 6 . The 
values of v are (a) 10~ 5 , (b) 1CT 5 , (c) 1(T 6 , (d) 10~ 7 , (e) 10~ 4 , and (f) 1CT 7 ; and the magnetic-diffusion time t v = ll/r) ~ 
2.52 x 10 10 , 2.52 x 10 8 , 2.52 x 10 8 , 2.52 x 10 4 and 2.52 x 10 4 , respectively. Dynamo action occurs in (a)-(d) but not in (e) 
and (f). 



it fluctuates significantly for a long time; the difficulty 
of pinpointing the dynamo boundary is a consequence of 
these fluctuations. 

The time series for the dynamo order parameter are 
obtained from those for E u and J5j; representative plots 
for these are shown, via red and blue-dashed curves, in 
Figs. 13 a), (b), (c), (d), (e), and (f) for a very large 
range of magnetic Prandtl numbers, namely, PrM = 
10 2 , 1, 10 _1 , 1CT 2 , 1CT 3 , and 1CT 6 , respectively. The 
values of v are (a) 10" 5 , (b) 1CT 5 , (c) 1CT 6 , (dj 10~ 7 , 
(c) 10 -4 , and (f) 10~ 7 ; and the corresponding values of 
the diffusion time scale t v = £q/t) ~ 2.52 x 10 10 , 2.52 x 
10 8 , 2.52 x 10 s , 2.52 x 10 4 and 2.52 x 10 4 , respectively. 
Clearly dynamo action occurs in Figs. [2ja)-(d) but not 
Figs.^e) and (f). By obtaining many such plots we can 
identify the dynamo boundary in the (Pr^ 1 , Pcm) plane 
as we discuss later. 

In the dynamo regime the shell- model kinetic and mag- 
netic energy spectra defined, respectively, by E u {k n ) = 
\u„\ 2 /k n and Eb(k n ) = \b n \ 2 /k n evolve as shown in 
Fig. [3l In particular, in Figs. [3ja), (b), and (c) for 
PrM = 10~ 2 , 1, and 10 2 , respectively, we show the evolu- 
tion of Eb(k n ) with time: the curves with red stars, green 
diamonds, blue hexagons, cyan circles, and magenta tri- 
angles, are obtained, respectively, for t = 1, 5, 10, 15, 
and 100; the analogs of these plots for E u (k n ) are given 
in Figs. G3d), (e), and (f). Note that the initial growth 
of Eb(k n ) occurs principally at large values of k n if PrM 
is large, i.e., we have a small-scale dynamo; this growth 
of Eb(k n ) moves to low values of k n as Pvm decreases; 
earlier studies [ill [2?| have observed similar trends but 
not over the large range of PrM we cover. As Ef,(k n ) 



grows, the velocity spectra are also affected but much less 
than their magnetic counterparts as can be seen by com- 
paring Figs. C2d), (e), and (f) with Figs. [3ja), (b), and 
(c), respectively. In all these plots the curves with black 
squares indicate E u (k n ) from the initial steady state for 
the GOY shell model; and the black lines with no sym- 
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bols show the K41 k n spectrum for comparison. From 
this line we see that the NESS that is obtained, once 
dynamo action has occurred, is such that both velocity 
and magnetic-field energy spectra display a substantial 
inertial range with K41 scaling; these inertial ranges are 
not large enough, at least near the dynamo boundary in 
our runs, for a reliable estimation of multiscaling correc- 
tions to the —5/3 exponent. If PrM — 1 then the scaling 
ranges in velocity and magnetic-field spectra are compa- 
rable; as PrM decreases (increases), the scaling range for 
the magnetic spectrum decreases (increases) relative to 
its counterpart in the velocity spectrum; these trends are 
clearly visible in the representative plots in Fig. [3J 

We return now to the identification of the dynamo 
boundary. A close scrutiny of the plots in Fig. [2] shows 
that the initial growth of P& is not monotonic. It is im- 
portant, therefore, to set a threshold value of the mag- 
netic energy E£: For a given pair of values for PrM and 
PeM, if Eb(t) > E^ for t > r c , where t c is the time at 
which the threshold value is crossed, we conclude that dy- 
namo action occurs; if not, then there is no dynamo for- 
mation. By examining the growth of Ef,(t) we can, there- 
fore, map out the dynamo boundary in the (Pr M ,Rey[) 
plane. The crossing time r c depends on PrM and PeM- 
Note that, if T c (PrM,PeM) > ^max, the length of time 
for which we integrate Eqs.Q and (|U), we would con- 
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FIG. 3: (Color online) Log-log (base 10) plots showing the time evolution of the magnetic-energy spectrum Et,{k n ) for rep- 
resentative parameter values at which dynamo action occurs: (a) Pryi = 1CP 2 , v — 10~ 7 ; (b) Pvm = 1, v = 10~ 5 , (c) 
Pr M = 10 2 , v = 10~ 5 ; analo gous plots for kinetic-energy spectra are shown in (d), (e), and (f), respectively. The curves with 
red stars, green diamonds, blue hexagons, cyan circles, and magenta triangles, are obtained, respectively, at t = 1, 5, 10, 15, 
and 100; the dissipation scale l d ~ 7.422 x 10" 4 in (b), (c), (e), and (f); l d ~ 4.779 x 10~ 4 in (a) and (d). Curves with black 
squares indicate velocity spectra before the seed magnetic field is introduced; the full black line shows a k~ 5 ^ 3 spectrum for 
comparison. 



elude, incorrectly, that no dynamo action occurs for this 
of values of PrM and PeM- In other words the dynamo 
boundary depends on i max ; we have checked this explic- 
itly in several cases. 

An important questions arises now: Is there a well- 
defined dynamo boundary in the (Pr^,Reyj) plane as 
i max — > oo? Earlier studies @, Us H3| have began to 
answer this question. They find that, if i max — t v , then 
a well-defined dynamo boundary is obtained. However, 
since they work with the MHD equations the error bars 
on this boundary are large and the range of values of 
PrM and PeM rather limited. 

The simplicity of our model allows us to carry out a 
systematic study of the dynamo boundary. We find that, 
at least in our shell model for MHD, we can obtain an 
asymptotic dynamo boundary (see Fig. 0]) if we choose 
P£ = Q.9E U , i.e., we conclude that dynamo action has 
occurred if Et(t) exceeds 0.9E U ; furthermore, if Eb(t) 
falls below 10 -35 we say that dynamo action will never be 
achieved. We continue the temporal evolution of Eqs.© 



and ([3]) till one of these criteria is satisfied. For all values 
of PrM and PeM that we have used we find that this i max , 
the run time required to decide whether or not dynamo 
action occurs, is several orders of magnitude lower than 
t^. We have also checked for several representative pairs 
of values for PrM and PeM that runs of length i max ~ r n 
do not change our conclusions about such dynamo action. 

The dynamo boundary that we obtain is shown in the 
stability diagram of Fig. [4j Red circles indicate param- 
eter values at which we obtain dynamo action whereas 
green stars are used for values at which no dynamo oc- 
curs. The most important result that follows from this 
stability diagram is that the boundary between dynamo 
and no-dynamo regimes is very complicated. It seems to 
be of fractal-type, with an intricate pattern of fine, dy- 
namo regions interleaved with no-dynamo regimes. This 
is especially apparent in the inset of Fig. [H which shows 
a detailed view of the stability diagram in the vicinity 
of the dynamo boundary. Earlier studies seem to have 
missed this fractal-type of boundary because they have 
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FIG. 4: (Color online) The dynamo stability diagram in the 
(Prj^ 1 , Reyi) plane: red circles indicate dynamo action; green 
stars are used if no dynamo occurs. The boundary between 
the two regions shows an intricate, interleaved pattern of 
fine, dynamo and no-dynamo regimes (see inset for a detailed 
view). We have drawn two black, dashed lines; the region 
above the upper one of these lines is predominantly in the 
dynamo regime; the area below the lower one of these lines is 
principally in the no-dynamo regime. 



not been able to examine the transition in as much detail 
as we have for our shell model. However, fractal-type 
boundaries between different dynamical regimes have 
been suggested in other extended dynamical systems; re- 
cent examples include the transition to turbulence in pipe 
flow [35l | and different forms of spiral-wave dynamics in 
mathematical models for cardiac tissue [3f|. In Fig. [4] 
we have drawn two black, dashed lines; the region above 
the upper one of these lines is predominantly in the dy- 
namo regime; the area below the lower one of these lines 
is predominantly in the no-dynamo regime. These two 
lines give an approximate indication of the error bars we 
might expect in the determination of the dynamo bound- 
ary in a study that cannot scan through points in the 
{Pr^,Rey\) plane as finely as we have. 

From Fig. Q] we see that the order parameter Ef,/E u 
jumps from a very small value in the no dynamo region 
to a value ~ 1 in the dynamo state. It is natural, there- 
fore, to think of the dynamo boundary as a nonequilib- 
rium, first-order boundary. In an equilibrium, first-order 
transition the order parameter shows hysteretic behavior 
if we scan through a first-order boundary by, say, chang- 
ing, at a finite rate, the field that is conjugate to the order 
parameter [30] . It is natural to ask if we see such hys- 
teretic behavior at the dynamo boundary. Indeed, we do, 
as we show in Fig. [5] where we cross the dynamo bound- 
ary by changing the amplitude fo of the forcing term in 
Eq.([3]). Figure [5] shows representative plots of the dy- 
namo order parameter E\,jE u versus fo; these illustrate 
the hysteretic behavior that occurs when fo is cycled at 
a finite, nonzero rate across the dynamo boundary; here 



Ptm = 10 -4 and v = 10~ 5 . As fo increases, E\,/E u fol- 
lows the blue, full line: it increases and then saturates; 
fluctuations are superimposed on these mean trends. If 
we now decrease fo, then Ef,/E u follows the red dotted 
line, and not the blue one, i.e., we have a hysteresis loop. 
The faster the rate at which we change fo the wider is 
the hysteresis loop as is known from studies of hystere- 
sis in spin systems [3(J. Here we increase fo in steps of 
1.0 x 10~ 3 from an initial value of 1.0 x 10~ 3 ; we keep 
fo constant for a time period 10 in Fig. OJa) and 1 in 
Fig. [Hb); the red, dotted-line segments of the hysteresis 
loops are obtained by decreasing f at the same rates as 
for the blue, full-line segments; the loop in the former 
case is narrower than in the latter. 

Given the analogy with first-order transitions that we 
have outlined above, it is natural to ask if nucleation-type 
phenomena [37| are also associated with dynamo forma- 
tion. It would be interesting to check this in a DNS of 
the MHD equations. At the level of our shell model, the 
best we can do is to try to see if, for a given Ptm, when 
we obtain a dynamo, the time required for dynamo ac- 
tion t c diverges as we approach the dynamo boundary. 
Our data are consistent with an increase of r c as we ap- 
proach this boundary from the dynamo side as shown by 
the representative plots in Fig. [51 However, it is hard 
to fit a precise form to the behavior of t c near the dy- 
namo boundary partly because of the complicated nature 
of this boundary which makes it difficult to estimate the 
position i?eMb reliably (the plot in Fig.[5Jis motivated by 
the form of Eq. (27) in Ref. " 



IV. CONCLUSIONS 

We have presented a detailed study of dynamo action 
in a shell model of turbulence [HI, 0, (25|. Our study 
has been designed to explore the nature of the bound- 
ary between dynamo and no-dynamo regimes in the 
(Prj^ 1 , Reyi) plane over a much wider range of Pryi than 
has been attempted in earlier numerical studies. The dy- 
namo boundary emerges as a first-order nonequilibrium 
phase boundary between one turbulent, nonequilibrium 
statistical steady state (NESS) and another [3g|. This 
point of view is implicit in earlier work, e.g., in studies 
of the Kazantsev dynamo (39j or in studies that view dy- 
namo generation as a subcritical bifurcation [3J HO, H]j . 
One of these studies [3l[ has remarked that when dynamo 
action "... is obtained in a fully turbulent system, where 
fluctuations are of the same order of magnitude as the 
mean flow ... the traditional concept of amplitude equa- 
tion may be ill-defined and one may have to generalize 
the notion of "subcritical transition" for turbulent flows 
We believe that the natural generalization is the 
nonequilibrium, first-order transition we suggest above. 
We have explored the explicit consequences of such a 
view in far greater detail than has been attempted hith- 
erto. In particular, the ratio Ef,/E u is a convenient or- 
der parameter for this nonequilibrium phase transition; it 
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FIG. 5: (Color online) Plots of the dynamo order parameter Eb/E u versus the forcing amplitude /o illustrating hysteretic 
behavior as /o is cycled across the dynamo boundary; here Prm = 10~ 4 and v — 10~ 5 . As /o increases, Et/E u follows the 
blue, full line; if we now decrease /o, then Eb/E u follows the red dotted line, and not the blue one, i.e., we have a hysteresis 
loop. We increase /o in steps of 1.0 x 10 -3 from an initial value of 1.0 x 10 -3 ; we keep /o constant for a time duration 10 in 
(a) and 1 in (b); the red, dotted- line segments of the hysteresis loops are obtained by decreasing /o at the same rates as for the 
blue, full-line segments. 
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FIG. 6: (Color online) Representative plots of lnr c versus Rem — Reya, for parameter values at which dynamo action occurs; 
here t c is in units of the time step 8t and fleMb is the estimated position of the dynamo boundary: (a) Pvm = 1 and (b) 
PfM — 5 x 10~ 4 . Note that the time t c required for dynamo action increases rapidly as we approach the dynamo boundary 
(the plot here is motivated by Eq. (27) in Ref. [13] )• 



shows hysteresis across the dynamo boundary like order 
parameters at any first-order transition; and nucleation- 
type phenomena also seem to be associated with dynamo 
formation. Last, and perhaps most interesting, we find 
that the dynamo boundary seems to have a fractal char- 
acter; this provides a natural explanation for the large 
error bars in earlier attempts to determine this bound- 
ary 0, EE Ull ■ Furthermore, this fractal-type boundary 
might well be the root cause of magnetic-field reversals 
discussed, e.g., in Refs. [29l. |42|. 

It is important to check, of course, that our shell-model 
results carry over to the MHD equations. This requires 
large-scale DNS that might well be beyond present-day 
computing capabilities if we want to explore issues like 
the possible fractal nature of the dynamo boundary. 
However, analogs of the hysteretic behavior we mention 



above have been obtained in DNS studies of the MHD 
equations [U 0, El[; hysteresis has also been seen in 
a numerical simulation that includes turbulent convec- 
tion [43|]. In some of these studies hysteretic behavior 
is obtained by changing the viscosity of the magnetic 
Prandtl number. We have obtained hysteresis by chang- 
ing the forcing; this change of forcing might be easier to 
effect in experiments than a change of the viscosity or 
magnetic diffusivity. 

To the best of our knowledge, earlier studies have not 
noted the increase in the dynamo-formation time r c as 
the dynamo boundary is approached from the dynamo 
side. We have suggested that this is akin to the increase 
in the time required to form a critical nucleus as we ap- 
proach a first-order boundary (3?| ■ It would be interest- 
ing to see if such an increase of r c can be obtained in DNS 



8 



studies of dynamo formation with the MHD equations. 
It is worth noting here that some DNS studies [Ty] have 
suggested that simulation times comparable to the dif- 
fusion time scale t v are required to confirm dynamo for- 
mation; by contrast our shell-model study yields dynamo 
action on a much shorter time r c , which increases as we 
approach the dynamo boundary. Perhaps the large simu- 
lation times required for dynamo action in full MHD sim- 
ulation might have arisen because these simulations have 
been carried out in the vicinity of the dynamo boundary. 

To settle completely whether the dynamo boundary is 
of fractal-type, very long simulations might be required 
to make sure that the apparent fractal nature is not an ar- 
tifact of long-lived metastable states. To make sure that 
our calculations do not suffer from such an artifact, we 
have carried out very long runs for representative points 
in the region of the dynamo boundary in Fig. [4j we have 
found that these long runs do not change our results. 
Furthermore, it is useful to check whether, instead of 
one dynamo boundary, there is a sequence of transitions, 
with more and more complicated temporal behaviors for 
the order parameter, as has been seen in the turbulence- 
induced melting of a nonequilibrium vortex crystal [44| . 
We have not found any conclusive evidence for this but, 
in the vicinity of the dynamo boundary, the order param- 
eter can oscillate for fairly long times (see, e.g., the green 
full curve in Fig.[T|). To decide conclusively whether these 
oscillations characterize a new nonequilibrium oscillating 
state, different from the simple dynamo and no-dynamo 
NESSs we have mentioned, requires extensive numerical 
studies that lie beyond the scope of this paper. 



In equilibrium statistical mechanics different ensembles 
are equivalent; in particular, we may determine a first- 
order phase boundary by using either the canonical or the 
grand-canonical ensemble. However, such an equivalence 
of ensembles does not apply to transitions between dif- 
ferent nonequilibrium statistical steady states (NESSs); 
examples may be found in driven diffusive systems [45| 
or in the turbulence-induced melting of a nonequilibrium 
vortex crystal j44|. Given that the dynamo boundary 
separates two turbulent NESSs, we might expect that 
this boundary might depend on precisely how the system 
is forced. Evidence for this exists already: For example, 
the dynamo boundary depends on whether a stochastic 
external force is used [7| or whether a Taylor-Green force 
is used [I6j ; furthermore, this boundary is different if the 
fluid is helical [21], as in most astrophysical dynamos. 

We hope our study of dynamo formation in a shell 
model for MHD will stimulate both DNS and experimen- 
tal studies designed to explore the first-order nature of 
the dynamo transition. 
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